- Plot vector geospatial data (points/polygons)
- Plot raster geospatial data (image)
- Extracting data from raster
- Making maps (base maps, compass rose, overlays, etc)
2020-02-24
We are researchers who want to answer the question "What range of environmental conditions do zebra tend to inhabit?"
General workflow:
Variants of this question/workflow are applicable to other study systems.
sf package for vector data
sp, rdgal for read/write, and rgeos for geometric operationstidyverseraster package for raster operations and more
sf objectssf (simple feature) objects are extended data.frames or tibblestibblesfc (simple feature columns) with useful components likebboxCRS (coordinate reference system) attributessf and tidyversesf functions begin with st_summary, plotsf methods for filter, arrange, distinct, group_by, ungroup, mutatemethods(class = "sf") ## [1] $<- [ ## [3] [[<- aggregate ## [5] anti_join arrange ## [7] as.data.frame cbind ## [9] coerce dbDataType ## [11] dbWriteTable df_spatial ## [13] distinct extent ## [15] extract filter ## [17] full_join gather ## [19] group_by group_map ## [21] group_split identify ## [23] initialize inner_join ## [25] left_join mask ## [27] merge mutate ## [29] nest plot ## [31] print raster ## [33] rasterize rbind ## [35] rename right_join ## [37] sample_frac sample_n ## [39] select semi_join ## [41] separate separate_rows ## [43] show slice ## [45] slotsFromS3 spread ## [47] st_agr st_agr<- ## [49] st_area st_as_sf ## [51] st_bbox st_boundary ## [53] st_buffer st_cast ## [55] st_centroid st_collection_extract ## [57] st_convex_hull st_coordinates ## [59] st_crop st_crs ## [61] st_crs<- st_difference ## [63] st_force_polygon_cw st_geometry ## [65] st_geometry<- st_interpolate_aw ## [67] st_intersection st_intersects ## [69] st_is st_is_polygon_cw ## [71] st_join st_line_merge ## [73] st_linesubstring st_make_valid ## [75] st_minimum_bounding_circle st_nearest_points ## [77] st_node st_normalize ## [79] st_point_on_surface st_polygonize ## [81] st_precision st_segmentize ## [83] st_set_precision st_simplify ## [85] st_snap st_snap_to_grid ## [87] st_split st_subdivide ## [89] st_sym_difference st_transform ## [91] st_transform_proj st_triangulate ## [93] st_union st_voronoi ## [95] st_wrap_dateline st_write ## [97] st_zm summarise ## [99] transmute ungroup ## [101] unite unnest ## see '?methods' for accessing help and source code
roads <- st_read(here("data", "enproads.shp"))
## Reading layer `enproads' from data source `C:\Users\mbuon\Desktop\QERM597-making-maps\data\enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID): NA
## proj4string: NA
class(roads)
## [1] "sf" "data.frame"
head(roads)
## Simple feature collection with 6 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.46841 ymin: -19.41752 xmax: 14.58964 ymax: -19.33077
## epsg (SRID): NA
## proj4string: NA
## TYPE LENGTH KM geometry
## 1 Track 0.7 1 LINESTRING (14.50502 -19.35...
## 2 Track 3.6 4 LINESTRING (14.49077 -19.36...
## 3 Graded 2.3 2 LINESTRING (14.4864 -19.417...
## 4 Track 7.1 7 LINESTRING (14.58498 -19.34...
## 5 Track 1.2 1 LINESTRING (14.58498 -19.34...
## 6 Graded 0.1 0 LINESTRING (14.52596 -19.33...
# Highlight parts
# 1. list column (aka sfc)
# 2. feature geometry (sfg)
# 3. feature (row)
# Default methods for objects plot(roads)
plot(filter(roads[,3], KM > 20))
rnaturalearthSpatial* objects using st_as_sfworld <- ne_countries(scale = "medium", returnclass = "sf") st_geometry(world) %>% plot()
Default data is Lat/Long
Namibia <- ne_countries(country = "namibia") Namibia_sf <- st_as_sf(Namibia) plot(Namibia_sf)
* Exercise: try for another country
#Uncomment here if using rnatural earth for the first time:
# devtools::install_github("ropenscilabs/rnaturalearthhires", force = TRUE)
# install.packages("rnaturalearthhires", repos = "http://packages.ropensci.org")
namibia_political <- ne_states(country = "namibia")
namibia_political_sf <- st_as_sf(namibia_political)
plot(namibia_political_sf)
st_crs("+proj=longlat +datum=WGS84") # Proj.4 string"
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(3857) # ESPG = 3857
## Coordinate Reference System:
## EPSG: 3857
## proj4string: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(3857)$units # check units
## [1] "m"
st_crs(4326)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(4326)$units
## NULL
st_crs(NA) # unknown (assumed planar/Cartesian)
## Coordinate Reference System: NA
Change the units!
units::set_units(a1, km^2) ## 824886.6 [km^2] units::set_units(a1, ft^2) ## 8.879005e+12 [ft^2]
zebra_data <- read.csv(here("data", "zebra.csv"))
head(zebra_data)
## X.1 X Name Date Longitude Latitude Speed Direction Temp
## 1 50 78 AG059 4/20/2009 15:40 15.89605 -19.12138 0 0 0
## 2 34 60 AG059 4/20/2009 13:01 15.91854 -19.17761 0 0 0
## 3 33 59 AG059 4/20/2009 12:02 15.91857 -19.17766 0 0 0
## 4 62 90 AG059 4/20/2009 17:00 15.80543 -19.08002 0 0 0
## 5 109 159 AG059 4/21/2009 1:40 15.91887 -19.17758 0 270 0
## 6 107 155 AG059 4/21/2009 1:00 15.91882 -19.17760 0 270 0
## Altitude PDOP ID DatePOS diff day
## 1 1123 2 AG059 4/20/2009 15:40 19 4/20/2009
## 2 1113 3 AG059 4/20/2009 13:01 59 4/20/2009
## 3 1130 3 AG059 4/20/2009 12:02 721 4/20/2009
## 4 1139 2 AG059 4/20/2009 17:00 17 4/20/2009
## 5 1133 2 AG059 4/21/2009 1:40 40 4/21/2009
## 6 1126 4 AG059 4/21/2009 1:00 20 4/21/2009
zebra_sf <- read.csv(here("data", "Zebra.csv")) %>%
dplyr::select(ID = Name, 4:6) %>%
mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
st_as_sf(., coords = 3:4, crs = "+init=epsg:4326")
st_transform(zebra_sf, 32733) #convert to UTM
## Simple feature collection with 11490 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 546746.9 ymin: 7852220 xmax: 691330.5 ymax: 7912185
## epsg (SRID): 32733
## proj4string: +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs
## First 10 features:
## ID Date timestamp geometry
## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501)
## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266)
## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260)
## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124)
## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269)
## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266)
## 7 AG059 4/21/2009 5:01 2009-04-21 05:01:00 POINT (596604.9 7879271)
## 8 AG059 4/21/2009 4:41 2009-04-21 04:41:00 POINT (596610.9 7879268)
## 9 AG059 4/22/2009 4:00 2009-04-22 04:00:00 POINT (596603.6 7879269)
## 10 AG059 4/22/2009 2:01 2009-04-22 02:01:00 POINT (596609.4 7879272)
countries <- ne_countries(scale=110) p <- ggplot() + geom_polygon(data=countries, aes(x=long, y=lat, group=group), color="white", lwd = .25) ggplotly(p)
st_writeZebra.csv and plot the geometry.roads <- st_read(here("data", "enproads.shp"), crs = "+init=epsg:4326") %>% #4326= coordinate system based on Earth's center of mass
st_transform(32733) #32733 = spatial reference for Nambia
## Reading layer `enproads' from data source `C:\Users\mbuon\Desktop\QERM597-making-maps\data\enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
st_geometry(roads) %>% plot
zebra_sf <- read.csv(here("data", "Zebra.csv")) %>%
dplyr::select(ID = Name, 4:6) %>%
mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
st_as_sf(., coords = 3:4, crs = "+init=epsg:4326") %>% #longlat, converting foreign object to SF object
st_transform(32733) #convert to UTM
ggplot() +
geom_sf(data=roads) +
geom_sf(data=zebra_sf, aes(color=ID))
How close or far from different road types do zebra move?
unique(roads$TYPE)
## [1] Track Graded Gravel Tar
## Levels: Graded Gravel Tar Track
#Filter roads based off type:
large_roads <- filter(roads, TYPE %in% c("Tar", "Gravel"))
small_roads <- filter(roads, TYPE %in% c("Graded", "Track"))
ggplot() +
geom_sf(data=large_roads, size=1.5) +
geom_sf(data=small_roads, size=0.6) +
geom_sf(data=zebra_sf, aes(color=ID))
# Find the minimum distance (in meters) of each point to a large road large_dist<- st_distance(y=zebra_sf, x=large_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds zebra_sf$large_road_dist <- apply(large_dist, 2, min) # Find the minimum distance (in meters) of each point to a small road small_dist<- st_distance(y=zebra_sf, x=small_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds zebra_sf$small_road_dist <- apply(small_dist, 2, min) head(data.frame(zebra_sf)) ## ID Date timestamp geometry ## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501) ## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266) ## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260) ## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124) ## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269) ## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266) ## large_road_dist small_road_dist ## 1 7.845860 3479.8987 ## 2 163.187201 241.8660 ## 3 156.746071 245.8446 ## 4 9.115608 1605.2154 ## 5 151.769751 227.9759 ## 6 151.358102 231.6984
ggplot(zebra_sf) +
geom_histogram(aes(large_road_dist, fill="large roads"), alpha=0.5) +
geom_histogram(aes(small_road_dist, fill="small roads"), alpha=0.5) +
scale_fill_manual(labels=c("large roads", "small roads"), values=c("blue", "orange"))
## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : C:/Users/mbuon/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max) ## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : C:/Users/mbuon/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max) ## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : C:/Users/mbuon/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif ## names : ndvi_mean_utm ## values : -558.117, 5303.021 (min, max)
## class : RasterLayer ## dimensions : 474, 1274, 603876 (nrow, ncol, ncell) ## resolution : 231.6564, 231.6564 (x, y) ## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : memory ## names : ndvi_mean_utm ## values : -0.0558117, 0.5303021 (min, max)
## ID Date timestamp geometry ## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501) ## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266) ## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260) ## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124) ## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269) ## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266) ## large_road_dist small_road_dist NDVI_mean ## 1 7.845860 3479.8987 0.2409661 ## 2 163.187201 241.8660 0.2682169 ## 3 156.746071 245.8446 0.2682169 ## 4 9.115608 1605.2154 0.2380762 ## 5 151.769751 227.9759 0.2419201 ## 6 151.358102 231.6984 0.2419201
Focus on ggmap - As of mid-2018, requires registering with Google and obtaining an API key - Requires providing a valid credit card (yikes!), though charges are nonexistent or at least very minimal - See the ggmap GitHub page for more information about API keys
ggmapggmapTwo steps:
This is my personal API key, which you can use for the purposes of this class!
register_google(key = "AIzaSyBRkw_wzDKuWZDikdD86Kmp1Sa6PzuCKFc")
To add a Stamen basemap, first define the bounding box for the basemap you would like to download.
zebra_bbox <- c(14, -20, 17.5, -18)
names(zebra_bbox) <- c("left", "bottom", "right", "top")
terr_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="terrain") ggmap(terr_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
wat_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="watercolor") ggmap(wat_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
To add a Google basemap, first define the center coordinates for the basemap you would like to download.
zebra_center <- c(15.8, -19)
names(zebra_center) <- c("lon", "lat")
sat_basemap <- get_googlemap(center=zebra_center, zoom=8, size=c(640, 640), scale=2,
maptype="satellite")
ggmap(sat_basemap) +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
Notice that get_googlemap() returns a square basemap tile, which then sets the plotting limits of your map.
To manually set different plotting limits, pass in an empty "base layer" and set x and y limits using ggplot().
ggmap(sat_basemap, maprange=FALSE, base_layer=ggplot(aes(x=1, y=1), data=NULL)) +
xlim(14.3, 17.4) + ylim(-20, -18) + xlab("lon") + ylab("lat") +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
RgoogleMaps
ggplot2 (boo)tif file as a raster brick, then display red/green/blue valuesgeocode function to find the lat/lon coordinates for a known location
geocode("Space Needle")my_basemap <- get_googlemap(center="Space Needle", zoom=14, size=c(640, 640),
scale=2, maptype="satellite")
ggmap(my_basemap)
## OGR data source with driver: ESRI Shapefile ## Source: "C:\Users\mbuon\AppData\Local\Temp\RtmpSWQUh6", layer: "ne_10m_rivers_lake_centerlines" ## with 1455 features ## It has 34 fields ## Integer64 fields read as strings: rivernum ne_id ## Warning in rgdal::readOGR(destdir, file_name, encoding = "UTF-8", ## stringsAsFactors = FALSE, : Dropping null geometries: 1453 ## although coordinates are longitude/latitude, st_intersection assumes that they are planar ## Warning: attribute variables are assumed to be spatially constant ## throughout all geometries